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Abstract. The anisotropic surface free energy, Andreev surface free energy, and 
equilibrium crystal shape (ECS) z — z{x, y) are calculated numerically using a transfer 
matrix approach with the density matrix renormalization group (DMRG) method. The 
adopted surface model is a restricted solid-on-solid (RSOS) model with "sticky" steps, 
i.e., steps with a point-contact type attraction between them (p-RSOS model). By 
analyzing the results, we obtain a first-order shape transition on the ECS profile 
around the (111) facet; and on the curved surface near the (001) facet edge, we 
obtain shape exponents having values different from those of the universal Gruber- 
MuUins-Pokrovsky-Talapov (GMPT) class. In order to elucidate the origin of the 
non-universal shape exponents, we calculate the slope dependence of the mean step 
height of "step droplets" (bound states of steps) («-(p)) using the Monte Carlo method, 
where p = {dz/dx, dz/dy), and (•) represents the thermal average. Using the result of 
the \p\ dependence of (n(p)), we derive a |p|-expanded expression for the non-universal 
surface free energy /eff(p), which contains quadratic terms with respect to The 
first-order shape transition and the non-universal shape exponents obtained by the 
DMRG calculations are reproduced thermodynamically from the non-universal surface 
free energy /cff(p). 



PACS numbers: 68.35.Md, 05.70.Np, 05.50.-fq, 68.35.-p, 05.10.-a 
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Figure 1. Hierarchy of models depending on the length scale. Example of a vicinal 
Si(OOl) surface. 



1. Introduction 

The surface free energy is one of the most fundamental quantities in surface science. 
The equihbrium crystal shape (ECS) (Fig. 1) is the shape of a crystal particulate with 
minimum surface free energy under equilibrium conditions[l]-[13]. A polyhedral ECS 
is obtained using the Wulff theorem from a polar graph of the anisotropic surface free 
energy through the Wulff construction. An ECS with several facets and curved areas is 
obtained using the Landau- Andreev method[7, 8], and this ECS is found to be similar 
to the Andreev surface free energy, where the work associated with step formation is 
ehminated from the surface free energy per projected area[8]. Since the ECS reflects 
the anisotropy of the surface free energy, studying the ECS corresponds to studying the 
surface free energy itself. 

Recently, we applied the restricted solid-on-solid (RSOS) model coupled with the 
Ising system (RSOS-I model) [14]-[16] to investigate the interplay between surface steps 
and adsorbatcs on a vicinal surface. The RSOS model[17] (Fig. 1 and Fig. 2) is an SOS 
(or Kossel crystal) [18, 19] model in which differences in height between nearest-neighbor 
(nn) surface sites are restricted to {0, ±1}. Prom statistical mechanical calculations 
using the RSOS-I model, a first-order shape transition is found to occur on the ECS 
profile. In addition to the shape transition, the "shape exponent" on the ECS seems to 
have values different from the universal Gruber-MuUins-Pokrovsky-Talapov (GMPT) or 
one-dimensional (ID) free fermion values [20] -[32]. 

Concerning the universal behavior of the GMPT, the universal form of the free 
energy is understood in terms of the terrace-step-kink (TSK) picture [20]- [23] (Fig. 1) 
for a vicinal surface as follows: 

f{p) = f{0) + ^p + Bp^ + O{p'), (1) 

where p represents the step density, 7 represents the step tension, and B represents the 
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(a) (b) 

Figure 2. Vicinal surface represented by the restricted solid-on-solid (RSOS) model. 

The difference in surface height between nearest neighbor sites is restricted to {0, ±1}. 
(a) Perspective view of surface tilted towards the (100) direction. A: terrace. B: step. 
C: kinks. D: adatom. E: ad-hole. F: island. G: negative island, (b) Top view of surface 
tilted towards the (110) direction. The gray lines represent surface steps and the filled 
squares represent points where adjacent steps collide. 

step interaction coefficient. A many-body system of non-overlapping linear excitations 
embedded in two dimensions is known to have the form of free energy as expressed in Eq. 
(1). The TSK picture has been confirmed by exact calculations of the free energy using a 
microscopic body-centered cubic solid-on-solid (BCSOS) model[24, 25]. Experimentally, 
the universal behavior expressed by Eq. (1) has been confirmed by observations of the 
"shape exponent" on the ECS profile [33]-[48]. The ECS near a facet is expressed by 
\z{x, y) — z{xc, yc)\ — A\x — Xc\^^^ , where {xc, He) are the coordinates of the facet contour, 
the power on the right hand side is the shape exponent in the normal direction 9^, and 
A is the "amplitude". From Eq. (1), and following some thermodynamic calculations, 
the value of 6n is obtained as 3/2, which is a GMPT universal value. 

In addition, based on recent developments in the study of non-equilibrium bunched 
steps [49]- [53], the values of the exponents in the profile of a bunched step are related 
to the force range of the effective step-step interactions on the non-equilibrium vicinal 
surface. The values of the exponents are thought to be related to the shape exponent 6'„ 
on the ECS. Prom these perspectives, the importance of studying the shape exponent 
is increasing. 

In order to establish the non-GMPT shape exponent on the ECS, a p-expanded 
expression for surface free energy with non-GMPT terms should be derived. In our 
previous work[15], we demonstrated the appearance of a short-range step-step attraction 
mediated by adsorbates, and we introduced the step-droplet picture [15]. Due to the 
complexity of interplay between surface steps and adsorbates, however, we could not 
derive a p-expanded expression for the non-GMPT surface free energy. 

The aim of the present paper is to establish the non-GMPT shape exponent on 
the ECS. In other words, the purpose is to find a mechanism to obtain a p-expanded 
expression for surface free energy with non-GMPT terms in a system with a short-range 
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step-step attraction. 

To obtain clear resnlts. we present a simple model: the RSOS model with a point- 
contact type step-step attraction (p-RSOS model) [54, 55]. Physically, the step-step 
attraction represents the transient bond formed by the spatial overlap of orbitals between 
atoms at the collision point of the adjacent steps (Fig. 2(b)). 

The paper is organized as follows. In §2, we present the definition of the p-RSOS 
model and show statistical mechanical calculations on the ECS, the equilibrium facet 
shape, and the non-universal shape exponents. In §3, a study on step droplets near 
equilibrium is described. Calculation of the mean step height of the step droplets (n) 
is carried out using a Monte Carlo method. In §4, we derive a p-expanded expression 
for the vicinal surface free energy, which contains non-GMPT terms. In §5, using the 
non-GMPT vicinal surface free energy, we thermodynamically reproduce the results for 
the p-RSOS model obtained by the transfer matrix method in §2. In §6, we present a 
summary and discussion. Finally, a conclusion is given in §7. 

2. Statistical mechanical calculations using the p-RSOS model 

2.1. Model Hamiltonian 

Let us consider the surface height h{i,j) at a site on a square lattice to describe 
surface microscopic undulations (Fig. 2). In the RSOS model[17], the height differences 
between nearest-neighbor (nn) sites are restricted to values of {1,0, —1}. We consider 
a point-contact type microscopic step-step interaction and refer to this model as the 
p-RSOS model. The Hamiltonian for the p-RSOS model can then be written as 



where e is the microscopic ledge energy, eint is the microscopic step-step interaction 
energy, and 6{a,b) is Kroneckcr's delta. The summation with respect to {i,j) is 
performed all over sites on the square lattice. The RSOS restriction is required implicitly. 
In the case of eint < 0, the interaction among steps becomes attractive. 

For a vicinal surface, we add the terms of the Andreev field[8] t] — {r]x,rjy) to the 
Hamiltonian Eq. (2) as an external field. The model Hamiltonian given in Eq. (2) for 
the vicinal surface then becomes 



"^p-RSOS — 



J2e[\h{i + l,j) - hii,j)\ + \hii,j + 1) - h{i,j)\] 



+ J2^i-t[S{\h{t + l,j + l)-h{t,j)\,2) 



+smz+i,j-i)-hii,j)i2)], 



(2) 





(3) 
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Figure 3. (a) Schematic diagram of the transfer matrix, (b) Graphical representation 
of the resultant matrix for the partition function. 



The partition function Z for the p-RSOS model is given by 

Z = ^ e-/3«vicinal (4) 
{h{i,3)} 

where /3 = l/k-oT, is the Boltzmann constant, and T is the temperature. The 
Andreev surface free energy /(T7) is the thermodynamic potential calculated from the 
partition function Z using 

where Af is the number of lattice points on the square lattice. Practically, calculation 
of Eq. (5) is not an easy task, because the entropy associated with the vast variety of 
zigzag structures of a surface step and by the parallel movement of steps is difficult to 
estimate. 

Recently, the numerical renormalization group method for one- dimensional (ID) 
quantum spin systems has been further developed to become the density-matrix 
renormahzation group (DMRG) method[56, 57, 58]. Though the method is approximate 
and numerical, it successfully reproduces known exact results with high precision. One 
means of extending the DMRG method to classical systems is by mapping a two- 
dimensional (2D) classical system to a ID quantum spin system[59] by use of the 
transfer matrix[60] together with the Suzuki- Trotter formula[61]. Such a method was 
developed by Nishino et al. for an infinite lattice, and is called the product-wavc-function 
renormalization group (PWFRG) method [62]-[64]. For the calculations in the present 
paper, we also adopt the PWFRG method. 

In order to apply the PWFRG method to the p-RSOS system, we construct the 
transfer matrix T{ti, t2, - ■ ■ , tN', s[, S2, • • • , s^) (Fig. 3(a)) using a 19- vertex model[65, 59] 
(Fig. 4). The partition function Z (Eq. (4)) is then rewritten in terms of T as 

Z = TT[f{t,,h, ■■■M; s[, 4, ■ • ■ , ^'iv)""] (6) 
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Figure 4. Mapping the solid-on-solid model onto the vertex model, (a) Surface 
configuration of the SOS model, (b) Mapping the configuration onto the vertex model, 
(c) Vertex model. 











±1 










±1 
















































±1 


±1 





±1 





1 




exp[-p8] 


exp[-pe] 


exp[-p£] 


exp[-pe] 


exp[-p(2 


e + 


±1 


±1 







±1 







1 




1 







±1 
















±1 




1 










1 







=Fl 



exp[-ps] exp[-ps] exp[-2pe] exp[-2pe] exp[-P(2s -h 8;^^ )] 



Figure 5. Explicit surface configurations corresponding to the 19 vertices and their 
statistical weights, where /3 = l/ksT. 



where N is the number of hnked vertices and M is the length of the system in the vertical 
direction in Fig. 3(b). The statistical weight of each vertex is shown in Fig. 5. For a 
vicinal surface, the statistical weight is multiplied by exp[^{(i + s')r]x + {s + t')r]y}/2\ 
with s, t, s' and t' having values of {0, ±1}. Then, by use of the statistical weight 
denoted by V{s,t; s',t'), the transfer matrix is expressed as follows (Fig. 3(a)): 

f{ti,t2, • • • , ijv; si, 4> • • • > Sjv) = XI ^2; 4> ?2) 

{9l},{92},- 

■ ■ ■V{qN-i,tN;s'j^,t'N). (7) 

In the limit M,N — >■ oo, only the largest eigenvalue of the transfer matrix A(A'") 
contributes to the partition function. The Andreev surface free energy, therefore, is 
obtained from Eq. (5) as 

/3/(77) = - lim ^ lnA(7V)^. (8) 

M,N^oo iV IVI 

The transfer matrix is diagonalizcd efficiently using the PWFRG method. In the 
PWFRG calculation, the number of so-called "retained bases" m is set from 7 to 37. 
The number of iterations for the diagonalization process is set to 200 ~ 10'^. 

We also calculate the surface gradient p — {px,Py) = {dz/dx,dz/dy) using the 
PWFRG method, where the surface gradient is defined as the thermal average of the 
height differences as follows: 

Pxitl) = {h{m + l,n) - h{m,n))d/ax, 
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Figure 6. Perspective views of tlie reduced ECS around tlie (001) facet calculated by 
the transfer matrix method with the PWFRG algorithm. (/3e)~^ = k-QT/e = 0.3. (a) 
p-RSOS model (eint/e = -0.5). (b) The original RSOS model (eint = 0). 




(a) (b) 



Figure 7. Normal line {PP') and tangential line {QQ') around the (001) facet on the 
reduced ECS. (a) Perspective view, (b) Top view. 

Py{r}) — {h{m,n + 1) — h{m,n))d/ay. (9) 

Here, (•) represents the thermal average, and and ay represent the lattice constants 
in the x and y directions, respectively {d — a^ — ay — 1). By sweeping the field ^t], we 
obtain curves for vs. f3r)x or py vs. Pr)y[14:, 15, 32, 54, 55, 65]. 

2.2. First-order shape transition on the ECS profile 

In Fig. 6, we show perspective views of the calculated Andreev surface free energy /{t]) 
divided by ksT around the (001) surface for k-oT/e — 0.3. Prom the thermodynamics of 
the ECS, the Andreev surface free energy and the Andreev field are related to the ECS 
coordinates [7, 8] by 

hVx, Vy) = ^z(x, y), rjx^ -Xx, rjy = -Ay, 

,.._f!Z), (10) 

where A is the Lagrange multiplier relating to the volume of the particulate. Eq. (10) 

implies that the surface shape of f{r]x,f]y) is similar to the ECS. We now introduce a 
reduced ECS Z = Z{X,Y) such that Z — Pf{r}x,r}y) — Xf3z{x,y), X — Prjx — —X^x, 
and Y = fSrjy = —Xf3y. 

In the case of ei^t < (attractive step-step interaction), large {111} facets appear in 
addition to the (001) facet and the {101} facets, because the (111) surface is energetically 
stabilized by eint- For comparison, we show the reduced ECS for the original RSOS model 
(eint = 0) calculated by the PWFRG method in Pig. 6(b). 
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Figure 8. Profile of tlic reduced ECS along the line PP'. (f)o — 7r/4. eint/e = —0.5. 
Broken lines represent metastable lines, (a) From right to left, k^T/e = 0.35, 0.36, 
and 0.37. (b) fceT/e = 0.36. The edge of the (111) facet is denoted by X^. (c) Original 
RSOS model (eint = 0). fceT/e = 0.3. 




Figure 9. Px vs. X along the normal line PP' at <^o = 7i'/4. (a) eint/e = —0.5. From 
right to left, ksT/e = 0.36, 0.361, 0.364, and 0.37. The broken line indicates in the 
metastable state, (b) Original RSOS model. ei„t = 0. k^T/e = 0.3. 



Along the line PP' shown in Fig. 7, we display the temperature dependence of the 
profile of the reduced ECS in Fig. 8. As seen from the figure, the (001) facet ends at 
(Xc, Yc) and the curved region between the (001) and (111) facets represents a first-order 
shape transition at the edge of the (111) facet {Xq, Xg)[54, 55]. 

Let us define the temperature Tf^i as the highest temperature at which the 
first-order transition occurs at the (111) facet edge on the ECS profile. We obtain 
ksTf i/e = 0.3610 ± 0.0005 for eint/e = —0.5 and 0o = 7r/4. For temperatures 
T < Tf^i, the first-order transition occurs at Xq{T) = Yq{T) (Fig. 8(b)). The surface 
slope changes abruptly from po = 1 to pi at Xq{T) for < pi < 1. In Fig. 9, 
we show Px{X,X) calculated by the PWFRG method along the line PP'. In the 
case of ksT/e — 0.36, p^ changes abruptly from pi = 0.349 ± 0.002 to po = 1 at 
Xq = Yq = 2.084 (> Xe = Fc = 2.0808 ± 0.0002). The values of p^ in the metastable 
state for 0.349 < p^ < 0.501 are shown by the broken line in the figure. 

Let us define the temperature T/_2 as the highest temperature at which the first- 
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Figure 10. Equilibrium facet shape (EFS) for X > and F > 0. Filled circles: 

{Xc,Yc) values calculated by the PWFRG method for k^T/e = 0.3. Open squares: 
{Yc,Xc) values. Solid lines: EFS for the 2D square nn Ising model (Eq. (11)) for 
fceT/e = 0.3. Dash-dotted lines: EFS for the Ising model for fceT/e = 0.361 (Eq. 
(11)). Dashed lines: Y = -X + 5.0. Dotted lines: Y = -X + 4.1551 (Eq. (12)). (a) 

Cint/e = -0.5. (b) Cint = 0. 



order transition occurs at the (001) facet edge. Below the curved area on the ECS 

between the (001) and (111) facets vanishes, and the (001) facet directly contacts the 
(111) facet. We obtain k^T^^^I^ = 0.3585 ± 0.0007 for ei„t/e = -0.5 and 0o = 7r/4. 

2. 3. Equilibrium facet shape 

In Fig. 10, we show the equilibrium facet shape (EFS) (Xc,Yc) of the (001) facet. The 
filled circles represent the (Xc, Yc) values obtained from the contour of Z{X, Y) = 0[32] 
calculated by the PWFRG method. Based on symmetry, we plot the {Yc, X^) values 
as open squares. The solid lines and dash-dotted lines in Fig. 10 represent the exactly 
calculated 2D ECS for the interface of the 2D square nn Ising model[66, 67, 28] for 
k-s,T/e — 0.3 and 0.361, respectively. The exact expression for the 2D ECS for the nn 
square Ising model is given by the following equation [68] -[71] (Appendix A): 

cosh(X,) + cosh(i;) = "^^^ (11) 

smh(/ye)) 

where ^ — l/k-QT. The dashed and the dotted lines in Fig. 10(a) represent the 
intersection between the (001) and (111) surfaces on the ECS for k-^Tje — 0.3 and 
0.361, respectively. They are approximately calculated by the following equation: 

y = -X + /3(2e + eint). (12) 

The facet shape obtained by the PWFRG calculations for k-^T je = 0.3 agrees well with 
the 2D ECS for the 2D nn Ising model except for the part truncated by the intersection 
line between the (001) and (111) surfaces. 
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Figure 11. Surface gradients calculated by the PWFRG method, eint/e = — 0.5. </io = 
7r/4 Filled circles: fceT/e = 0.36. Filled triangles: k^T/e = 0.37. = = 2.0808 
for kBT/e = 0.36, and Xc = Yc = 2.0051 for teT/e = 0.37, (a) ln(p^) w. ln[r/V2]. 
Solid line: In(p^) = 0.981n[r/\/^] + 4.4. Broken line: ln{p^) = 0.500 ln[r/\/2] - 0.088. 
(b) ln|p^ -py\ vs. ln[|i|/\/2]. Solid line: ln\p^ - Py\ = 2.96 ln[|i|/\/2] +4.4 Broken 
line: In \p^ - Py\ = 1.8 ln[|i|/v^] - 0.79. 

The step tension (or the interface tension) 7(0) is given by[32] (Eq. (C.2)) 

^7(0) = COS + ^0 sin 0. (13) 

Therefore, the agreement between the EPS and the 2D ECS for the 2D Ising model also 
implies an agreement between the step quantities such as the step tension and the step 
stiffness (7(0) = 7(0) + 9^7(0)/50^) for the RSOS model and the interface quantities 
such as the interface tension and the interface stiffness for the 2D Ising model [66] 
(Appendix A). For T < Tf^2, the first-order shape transition occurs at the (001) facet 
edge. In this case, a "step" actually corresponds to a "giant step" with a height nd 
{d = 1). In the large step-height limit, lim„_>oo7n('?'"/4)/n converges to \/2(e -|- eint/2), 
which is smaller than 7(7r/4)ising. 

Let us now calculate the approximate value of For T < the EPS has the 
shape of the rounded square truncated by the intersection line between the (001) and 
(111) surfaces, which is expressed by Eq. (12). Por T — Tf^2, this line contacts the EPS 
at {XI, XI) as the tangent line. Then, from Eq. (11) and (12), we have 



By solving Eq. (14), we obtain e/fceTj."^ ^ 2.78778 or kBTj^e ^ 0.358709, which is 
consistent with the PWPRG calculated value of k-oT f 2/^ (§2-2). 

2.4- Non-universal shape exponents 

Let us assign r and 0o to the line PP' for X > so that X = X^ + r cos 0o and 
Y — Yc-\-r suKpQ, where 4>q = n/Ais the tilt angle of the line PP' at (X^,, Yc) relative to 




(14) 
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the crystal axes (Fig. 7 (b)). The "normal shape exponent" 6^ is defined as the shape 
exponent[32] along PF on the ECS profile such that \Z{X{r),Y{r)) - Z{Xc,Yc)\ = 
An{<po)f^^" (0 < r), where we refer to the coefficient An{<po) as the "normal amplitude". 
Similarly, the "tangential shape exponent" 6t and the "tangential amplitude" At{(f)o) 
along are defined such that \Z{X{t),Y{t)) - Z{Xc,Yc)\ = ^(00)1^1^*, where t is a 
parameter assigned to the fine QQ' as follows: X — Xc — t sin (j)^, and Y — Yc + t cos 0o- 
Recalling that 0o — t^I^i we express {pr,Pt) in terms of {px,Py) as follows: 

Pr = \/2px, Pt = {Py - Px)l^- (15) 

We show the logarithm of p^ and pt in Fig. 11(a) and (b), respectively. 

First, we study the shape exponents for k-oT/t = 0.36, which is an example of 
Tf,2 < T < Tf^i. By fitting the data in the range —7.8 < ln[r/\/2] < —6.5 in Fig. 
11(a) to the linear function Aq + Ailn[r/\/2] by the least squares method, we obtain 
Ai = 0.98 ± 0.03, and Aq = 4.3 ± 0.2. Similarly, by fitting the data in the range 
-4.75 < ln[|i|/V2] < -4 in Fig. 11(b) to A'q + A[ ln[|i|/v^], we obtain A[ = 2.96±0.08, 
and Aq — 4.4 ± 0.2. These values give the shape exponents and amphtudes as 
9n = 1.98 ± 0.03, 9t = 3.96 ± 0.08, A(vr/4) = 40.6 ± 0.2, and A(7r/4) = 7.2 ± 1.0. 

Both exponents disagree with the GMPT universal values of 6n = 3/2 and 
6t = 3 [32]. Consequently, we conclude that the profile near the (001) facet contour 
for Tj2 < T < Tj i and (pQ = it/A shows non-GMPT behavior in the limit Pr, Pt — ^ 0. 

Next, we study the shape exponents for k-QT/e = 0.37, which is an example of 
T > Tji. As seen from Fig. 11(a), the slope of Inp^. crosses over from the larger 
value to the smaller value as Inr decreases. In Fig. 11(a), the broken fine represents 
Inpx — 0.51n[r/-\/2] — 0.088, and this fine is determined by fitting to the three lowest 
points. Therefore, in the hmit r — >■ 0, we have 9n = 1.5 and ^„(7r/4) = 0.73 ± 0.04. For 
Inpt, we fitted the data in the range -5 < ln[|t|/\/2] < -4.3 to Aq + A[ ln[|t|/v^] and 
obtained A[ = 1.8 ± 0.3 and A'q = -0.79 ± 0.08. These values lead to Ot = 2.8 ± 0.3 and 
At{Ti/4:) = 0.062 ± 0.018. 

From the values of 6^ and 6t for k^T /e = 0.37, we conclude that the profile near 
the (001) facet contour for T > Tf^i and (pQ = 7r/4 behaves like the GMPT universal 
profile in the limit pr, Pt ^ 0. 

We obtain the GMPT amphtudes using the equations[32] 

2v/2^^ 1 
M^o)-^^, M^o)-^^^^y (16) 

where we use the universal relation [2 9] 

/3i?(0o)=7rV[6/37(0o)]. (17) 

In the temperature region T > from the results shown in §2.3, the step quantities of 
a single step is found to be well described by the interface quantities of the 2D nn Ising 
model. The exact expressions for the interface tension 7(0) and the interface stiffness 
7(0) in the 2D nn Ising model for the case of (po = tt/A are as follows[66, 67, 28] (Eq. 
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Figure 12. vs. X along the normal line PP'. (po — 7r/4. (a) ksT/e = 0.36. 
ei„t/e = -0.5. = Fc = 2.0808. Filled circles: PWFRG results. Thin solid line: 
= 0.4432V^ - for X > X^. Broken line: p,, = 80.4X- 167.3 (b) feT/e = 0.37. 
eint/e = -0.5. = Yc = 2.0051. Filled triangles: PWFRG results. Thin solid line: 
p^ = 0.4421VX - Xc for X > X^. Broken line: p^ = 0.9154^^ - Xc for X > Xc- 
(c) Original RSOS model, eint = 0. kBT/e = 0.36. Thick solid line: PWFRG results. 
Thin solid line: = 0.4432 - Xc for X > Xc- 



(A.4)): 



For kBT/e 



/37(^)lsing 



\/2 cosh 



coshVf) 
2sinh(/3e) 



^7(^)ising = v^tanh [/37(-)isi„g/^ 



(18) 



2.08075, and /37(7r/4)ising 
i 2.0051, and /37(7r/4)ising 



(/3e)-i = 0.36, we have /37(vr/4)isi„g/V2 ? 
1.371, and at ksT/e = 0.37, we have /57(7r/4)ising/\/2 
1.364. Using these values for the interface stiffness, we can obtain the GMPT amphtudes 
from Eq.(16). 

In Fig. 12, we show Px{r) along the normal line PP'. The solid curve in Fig. 12 
shows Px{r) calculated from Eq. (15) together with the GMPT amplitude determined 
using Eq. (16). As seen from the figure, the solid curve disagrees with the values 
calculated by the PWFRG method not only for k^T/e = 0.36 (Fig. 12(a)) but also for 
ksT/e = 0.37 (Fig. 12(b)). For ksT/e = 0.36, with p^. < 0.08, fitting the values obtained 
by the PWFRG method to = A + A^X gives Ao = 167.3 ± 0.3 and Ai = 80.4 ± 0.2. 
This is shown as the broken line in Fig. 12(a). For k^T /e = 0.37, the broken line in 
Fig. 12(b) represents Px — 0.9154:\^X — X^ based on the value of Aq obtained from the 
straight fine in Fig. 11(a). The curve agrees with p^ calculated by the PWFRG method 
for small X — Xc- 



2.5. Vicinal surface tilted towards the (100) direction 



Let us now consider the profile along the curve OP in Fig. 13(a). We show the calculated 
Andreev surface free energy and the Px-X curve in Fig. 13(b) and (c), respectively. For 
00 = 0, we obtain p^ and Py from Eq. (10) and Eq. (16) and the GMPT universal shape 
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Figure 13. Profile of tlie reduced ECS along the curve OP in the p-RSOS model, 
fint/e = —0.5. (a) Schematic diagram of the ECS and the curve OP (broken line), (b) 
f/kBT (Z) vs. X. From right to left, feT/e = 0.2, 0.3, and 0.4. (c) Surface slope 
vs. X. fceT/e = 0.2. Xc = 4.9865, Yc = 0. Broken line: calculated using Eq. (19) 
withEq. (20). 



exponents as 



Px = 



TT 



1 



Y=0, 



Py = 



\X=Xr. 



(19) 



7r/37(0) 

where Fc = and = 4.9865 ± 0.0005. The interface tension and stiffness in the Ising 
model are obtained exactly from Eq. (13) and Eq. (A. 4) as follows: 

cosh^(/3e) 



/37(0)lsmg 

/37(0)lsmg 



cosh" 



- 1 



sinh(/9e) 

sinh [/37(0)ising] . (20) 

The value of Xc agrees well with /37(0)ising- Px calculated using Eq. (19) with Eq. (20) 
is plotted as the broken line in Fig. 13(c), and it can be seen that it closely matches the 
Px-X curve calculated by the PWFRG method for small p^- 

Therefore, we conclude that the vicinal surface tilted towards the (100) direction 
shows typical GMPT universal behavior. 



3. Origin of non-universal behavior: Step droplets 

3.1. Vicinal surface free energy 

When the surface slope is chosen to be an external variable instead of the Andreev field 

r], the thermodynamic function of the surface becomes the vicinal surface free energy 
f{p). Using the relationship between the Andreev surface free energy and the vicinal 
surface free energy[8], we obtain f{p) from /(ry) as 

f{p)-f{ri) + VP- (21) 
Namely, using the notations on the reduced ECS, we have 

(3f{p)^Z + Xpx + Ypy. (22) 
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Figure 14. Surface free energy per projected area /(p) (vicinal surface free energy). 

(a) p-RSOS model, ejni/e = —0.5. Px = Py = P- Solid line: f(p.p) at k^T/e = 0.36. 
Dashed line: f{p,p) at kBT/e = 0.36 for the metastable state. Open circle: 
at k-QT/e = 0.36. Open squares: /(0,0) and /(1,1) at k-QT/e = 0.3. (b) The original 
RSOS model (ei„t = 0). Solid line: f{p,p) at feT/e = 0.3. 




Figure 15. Thermal step bunching and coexistence of two surfaces, (a) Schematic 
diagram of vicinal surface free energy f{p) = f{p,p). Filled circle: f{po) = ,/(!)■ 
Broken line: tangential line connecting and {po,f{Po)) with a slope of Rq 

(Eq. (24)) [15]. (b) Schematic diagram of profile of the reduced ECS. R = VX^ + Y^. 
Broken line: shape in the metastable state. The (111) facet edge, where the first-order 
transition occurs, is indicated hy Rq. A surface with a slope po coexists with a surface 
with a slope pi at Rq. (c) Step bunching when two surfaces coexist. Dotted line: mean 
surface slope (= p). Broken line: local surface with slope pi. 



In Fig. 14(a), /(p) in the p-RSOS model (eint/e = -0.5) calculated by Eq. (22) is 
shown, and f{p) in the original RSOS model is shown in Fig. 14(b). For k-QT /e = 0.3, 
only /(0,0) and /(1, 1) are obtained by the PWFRG calculations, because a vicinal 
surface with a regular train of steps in the region < p < 1 is thermodynamically 
unstable. For kBT/e — 0.36, a f{p) curve is obtained in the region < Px < 0.349, and 
a curve for the metastable state in the region 0.349 < Px < 0.501. For k^T /e — 0.3 and 
0.36, /(1, 1) is well approximated by /(1, 1) = 2e + eint- 
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Figure 16. Snapshots of the simulated vicinal surfaces. Top view, eint/e = —0.5. 

240-\/2 X 240-\/2. Surface height is represented by brightness with 10 gradations, 
with brighter regions being higher. The black areas next to the white areas represent 
the higher terraces than the ones by one because of the finite gradation, (a) Initial 
configuration. Number of steps = 24. (b) fcsT/e = 0.1. 5 x 10*^ MCS/site. (c) 
knT/e = 0.35 (number of steps = 10). 4 x 10* MCS/site. 



3.2. Thermal step bunching 

The first-order shape transition on the ECS profile around the (111) facet (Fig. 15(b)) 
leads to the coexistence of two surfaces in equilibrium (Fig. 15(a)) [15, 16, 54], and we 
refer to this process as thermal step bunching[14, 15]. Let us consider the free energy 
fip) {p = \p\) along PP' in Fig. 7 for a mixture of surfaces with slopes of Pq and pi so 
that the mean slope is p, as shown in Fig. 15(c). The free energy of the simple mixed 
system is described as 



f{p) = xif{pi)+xof{po) 

Xi+ Xo = I 



(23) 



where Xq and Xi represent the fractional areas of the surfaces with slopes of po and pi, 
respectively. Since xq — {p — Pi)/(po — Pi) and xi — (po — p)/ipo — Pi), f{p) can be 
rewritten as 



f{p) = /(Pi) + 



[/(Po) - fiPi)] 



(P-Pi)- 



(24) 



Po -Pi 

Eq. (24) describes the co-tangent line which contacts f{p) at points {po,f{Po)) and 
{pi, f{pi)) (Fig. 15(a)). For pi < p < po, therefore, the free energy of the mixed surfaces 
is lower than that for a homogeneous surface. Moreover, since R — /3df{p)/dpr, Eq. 
(24) is equivalent to Z(R,) = l3~f(R,) = I3f(p) - R,p = l3f(R)\R^R^,- = l3f(R)\R^R^,+, 
where p = 1/kBT (Fig. 15(b)). 



3.3. Mean step height 

3.3.1. Giant steps: T < T/_2 We consider that the origin of the non-GMPT shape 
exponents is the formation of local "step droplets" [14, 15], similar to the formation of 
clusters in a gaseous system near the transition temperature. In order to form a clear 
image of such step droplets, we investigate the step dynamics near equilibrium using 
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a Monte Carlo method for a vicinal surface tilted towards the (110) direction at low- 
temperature. 

Initially, steps numbering A'^stcp are spaced evenly on a surface with an area of 
240^2 X 240^2 (Fig. 16(a)). The left side of the image is higher than the right side 
by an amount equal to iVstep- Periodic boundary conditions are imposed in the vertical 
direction in Fig. 16. 

To study the time evolution of the step configuration, we adopt a simple Metropohs 
algorithm without any driving force to simulate crystal growth. We randomly choose 
a site and allow its height h{i,j) to increase or decrease with equal probability. 

Then, if the RSOS restriction is satisfied, the height is updated by the Metropolis 
algorithm with a probability P described by 



where AE{i,j) — E{h{i,j) ± 1) — E{h{i,j)). The energy E{h{i,j)) is calculated using 
the p-RSOS Hamiltonian shown in Eq. (2). 

In the present Monte Carlo simulation, only non-conserved attachments and 
detachments of atoms are taken into consideration. Other effects that occur on a 
real surface are ignored, such as surface diffusion[3], clectromigration[72]-[74], the 
Schwoebel effect[75, 76], impurity effects [77] -[79], strain effects[80]-[84], and the effect 
of surfactants [85]- [86]. 

For T < Fig. 16 shows the formation of giant steps, similar to 

"step faceting" [87, 4]. As mentioned in our previous paper [55], at sufficiently low 
temperatures, step zipping occurs (Fig. 16(b)), whereas step unzipping seldom takes 
place. The zipping process corresponds to successive sticking together of steps, starting 
from the colliding point of adjacent steps, similar to the action of a zip fastener. At 
slightly higher temperatures where unzipping occurs more frequently, the vicinal surface 
reaches an equilibrium configuration (Fig. 16(c)). All the steps join together to form a 
single giant step, whose edge has the appearance of a smooth (111) surface. 

3.3.2. Step droplets: < T < Tf^i In this temperature regime, the fiuctuations on 
the step meandering become larger. A snapshot of the surface ior U-qT /e — 0.36 is shown 
in Fig. 17. In contrast to the low-temperature surface structure, giant steps are seen to 
exist locally, and become larger and smaller dynamically. We refer to such local giant 
steps as "step droplets" . Based on the structure of giant steps at low temperature, we 
now use Monte Carlo simulations to determine the size of step droplets and the mean 
number of steps (n) contained in giant steps. 

In order to calculate (n), we divide the square lattice into two sub-lattices, which 
are represented by open circles and closed circles in Fig. 18, respectively. Sweeping 
i from left to right for fixed j on the A-sublattice, we define the beginning of a step 
droplet at (i, j) as the position where successive height changes begin to occur. If these 




(A£;(i,i)<0), 

(A£;(i,j) >o). 



(25) 
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Figure 17. Snapshot of simulated vicinal surfaces. Top view, eint/e = —0.5. 4 x 10^ 
MCS/cite. 240\/2 x 240\/2. Number of steps: 24. Surface height is represented by 
brightness with 10 gradations, with brighter regions being higher. ksT/e = 0.36. 



n{i,j)=2 

.O. . O.vO.vO . o. 




I Top View 
}l — O— O-i Side View 




(a) 



n{i,j)=3 




I Top View 
— O— O-i Side View 




(b) 



Figure 18. Schematic illustration of bunched step. Open circles: A-sublattice points. 

Closed circles: B-sublatticc points. Upper figure: top view. Bottom figure: side view 
of the j-th row (shaded area) in the top view. Solid line (top view): surface steps, (a) 
n{i,j) = 2. (b) n{i,j) = 3. 



successive height changes ends at (i + n,j) in the A-sublattice, the size of the step 
droplet is taken to be n{i,j). The mean size of the step droplets (n) is then defined by 

(n)-^ EEl^(^'^')l' (26) 

where = Nj — 240 is the number of lattice units in the i- and j- directions, 
respectively. iVdropiet is calculated using A^dropict = Ejii ^droplet where iVdropiet(j) 
represents the number of giant steps (step droplets) in the j-th row of the square lattice 
shown in Fig. 18. 
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Figure 19. Slope dependence of mean step-droplet size. Monte Carlo calculation. 
Lattice size: 240v^ x 240\/2. ei„t/e = -0.5. Open circles: fceT/e = 0.36, 8 x 10^ 

MCS/site. Open triangles: fceT/e = 0.37, 4 x 10® MCS/sitc. Open reversed triangles: 
knT/e = 0.36, a^t = (original RSOS model), 2 x 10® MCS/site. 

In Fig. 19(a), we plot (n) calculated using the Monte Carlo method against Pr, as 
expressed in Eq. (15). As can be seen, (n) increases linearly with pr. Fitting the data to 
(n) = Aq + AiPr using the least squares method, we obtain the values Aq = 1.00 ± 0.02 
and Ai = 2.26 ± 0.08. To confirm the linearity, ln[(n) — 1] vs. Inpr is plotted in Fig. 
19(b). Fitting the data for > 0.03 to ln[(n) - 1] = Aq + Ai In yields = 0.84±0.04 
and Ai — 1.00 ± 0.02. The upper fitted line in Fig. 19(b) is seen to agree with the data 
for pr < 0.03 (open circles). 

3.3.3. T > 1 In this temperature regime, the GMPT universal behavior is expected. 
As seen from Fig. 19, however, the GMPT universal behavior such as that for the 
original RSOS model (reversed open triangles) is not observed for /ceT/e = 0.37 (open 
triangles) . 

Fitting the data for the original RSOS model to (n) = Aq + AiPr + ^2Pr yields 
^0 = 1-001 ± 0.008, Ai = -0.07 ± 0.09, and A2 = 2.2 ± 1.0. These values suggest that 
Aq — I and Ai — for the GMPT universal system in the limit Pr — > 0. 

For the p-RSOS model, (n) seems to increase linearly with pr for k-oT /e — 0.37. 
Fitting the data for k-B.T/e = 0.37 in Fig. 19(a) to (n) = A'^ + A'^Pr yields A[, = 
0.996 ± 0.04 and A'-^ = 1.45 ± 0.03. Fitting the data for ksT/e = 0.37 and pr > 0.03 in 
Fig. 19(b) to ln((n) - 1) - A'^ + A'llnpr yields A'^ = 0.58 ± 0.5 and A'l = 1.09 ± 0.07. 

4. p-expanded expression for non-GMPT vicinal surface free energy 
4.1. Breakdown of homogeneous ID fermion system 

Let us consider the mean running direction of steps as the direction of time flow. If we 
consider a coarse-grained vicinal surface of the original RSOS model, the steps can be 



Non-universal equilibrium crystal shape 



19 



lo OO O OO O O I 

/2(P2) f<„>M 
100 0000 0000 001 lo 00 O O OOOOOO I 

w 

looo 000 oooood~| ^ 



Figure 20. Schematic illustration of ID many-body systems with boson n-mers, and 
their mixture, (see §4.2.) 



regarded as continuous lines describing the space-time trajectories of Brownian particles 
with hard-core repulsion. If the particle density is low, the Hamiltonian for the transfer 
matrix in the continuous system is known to be described by ID free fermions[21]- 
[23], [29], and the free energy is expressed by Eq. (1). 

In the case of eint 7^ 0, since two fermions can not be present at the same site at the 
same time, the transfer matrix for the p-RSOS model for a continuous system cannot 
be written as the Hamiltonian for interacting fermions. Using the transfer matrix of the 
RSOS model, which is similar to the p-RSOS model, den Nijs and Rommels mapped the 
RSOS model to a ID quantum spin system[59]. Since the domain walls in the ordered 
phase correspond to the steps on the surface, the Hamiltonian for the transfer matrix in 
the continuous system is expressed in terms of interacting impenetrable bosons[88, 89]. 

In the case of eint < 0, bound states such as molecules or clusters might be formed 
between nn bosons [90] depending on the strength of the attraction between the bosons. 
That is, 7 and B in Eq. (1) can vary depending on the properties of the bound state, 
such as the number of particles in a cluster. In contrast, in a ID free fermion system, 
since 7 and B are the microscopic quantities assigned to a fermion monomer in the 
continuous model, changes in these parameters are not allowed. In this way, the ID 
fermion picture breaks down, and we must take inhomogeneous effects such as cluster 
formation in the Bose gas into consideration (a similar argument is given in §3.2). 

4.2. Expression for non-GMPT vicinal surface free energy 

We regard a giant step with a size n as a boson n-mer formed by microscopic step-step 
attraction eint (< 0). As time passes, these n-mers may grow or shrink in a similar 
manner to crystal clusters in the gas phase. Recalling the results described in §3, we 
assume a system of ID bosons made up of a mixture of different boson n-mers as shown 
in Fig. 20. We refer to this inhomogeneous picture of a vicinal surface as the step-droplet 
picture. 
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Let us now consider such a many-body system of boson n-mers (Fig. 20). We 
denote the density of giant steps with height n by p„. Wc assume that the interactions 
among n-mers are short range. The free energy fn{Pn) of a vicinal surface made up of 
such giant steps then becomes the GMPT universal type since 

UPu) = /(O) + 7n(0)Pn + B^{4^)pI + CMp'u + 0{pI), (27) 

where 7^(0) represents the step tension of a giant step, Bn{(f)) represents the step 
interaction coefficient among giant steps, C„(0) (< 0) represents the three-body 
interaction coefficient among giant steps[91], and 0, similar to the case for a single 
step (Fig. 7), represents the angle between the mean running direction of the giant step 
and the y-axis. In the limit p„ — )■ 0, converges to 0o. 

Recalling that pn = Ns,n/ L = \p\/dn = \p\/{ndi), where Ns^n is the total number 
of giant steps with size n, L is the length of the projected area of the vicinal surface, 
and dn is the height of the giant steps, we obtain a |p|-expanded expression for fn{Pn) 
as fn{p) (Appendix B, Eq. (B.l)). We then expand 7„(0) and Bn{(t>) with respect to n 
around n — 1. After some calculations (Appendix B), we obtain 



U\P\) = /(o) + 



7i 



-n 



\P\_ 
di 



7i(0)-7i'^(0) + 
5i(0) \P\' 



7f^(0) 



(2), 



dl 



l + (n 



' 71 (0) 



1, 
+ 2(- 



1)2#(0) 



7i(0) 



-1 



a(0) \p\ 



+ 0{\p\'), 



d\ 

where 7^'"'' (0) and 7^™''' (0) are defined by 



(m) 

7i 



(0) 



In 



(& /n 



#)(0) 



a-(7„(0)/n) 



(28) 



(29) 



and 7n(0) = 7n(0) + <9^7n(0)/f^0^ represents the stiffness of a giant step. 

Since n is difficult to estimate, we approximate Eq. (28) by replacing n with (n(0)). 
Recalling the results for {n{(j))) obtained by the Monte Carlo calculations in §3.3, we 
expand {n{4>)) with respect to \p\ around |p| = 0. Namely, 



n, 



(m) 



(n(0)) 



l + 4'^(0)|p| + -n, 



d\p\ 



(30) 



IPI=o+ 



Finally, substituting Eq. (30) into Eq. (28), and after some calculations, the non- 
GMPT vicinal surface free energy can finally be expressed as 



/eff(p) = f{n){p) = /(O) + 7l(0)|^ + ^eff(0)|pr 
+Seff(0)|p|^ + Cefr(0)|p|^ + C»(p^), 



(31) 
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where 



+ 



(32) 

(33) 



CM) 



f (0)7f)(0) + 3rz«(0)4^)(0)7f^(0)] 

Cn{4>) 



4+f (^) 



+ 



(34) 



dl \ 7i(0) I d\ ■ 
In Eq. (31), the term |pp appears. In addition, Bes{4>) can be larger or smaller than 

The microscopic quantities e and e^^t do not appear expUcitly in Eq. (31). They do, 
however, seem to determine mesoscopic quantities such as step tension, stiffness, and 
mean height. The mean step height depends on the character of the step-step attraction 
(see Fig. 19(a)). In fact, the slope dependence of the mean step height for the RSOS-I 
model seems to be different to that for the p-RSOS model. A detailed study on the 
RSOS-I model will be reported elsewhere [92]. 

4-3. Thermodynamic expression for non-GMPT shape exponents 

In order to derive the non-GMPT shape exponents thermodynamically, we write the 
|p I -expanded form of the vicinal surface free energy as 

hip) = /(O) + 7(0)|P| + B^{4>)\p\^ + 0(|p|^+^), (35) 
where C > 1- We can then express the shape exponents thermodynamically as (Appendix 
C) 



^, = 2C/(C-i), A(0) 



(kBT) 



C+i 



c-i 



for the shape exponent along the tangential line {QQ' in Fig. 7), and 



en = (/{(- 1), M<P) 



On 



CBd<l>) 



c-i 



(36) 



(37) 



for the shape exponent along the normal line {PP' in Fig. 7). From Eq. (36) and Eq. 
(37), we obtain the following relationships: 

0t = (38) 



A(0) 

,A(0) 



kBT 
27(0) 



(39) 



5. Application to p-RSOS model 

5.1. Tf,2<T<Tf,^ 

We compare the results of Eqs. (36) and (37) for the non-GMPT free energy Eq. (31) 
with the numerical results in §2 obtained by statistical mechanical calculations. 
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Figure 21. Schematic diagram of slope dependence of r. Solid lines: r{pr) of the 
ECS. Broken lines: Formally obtained r{pr) using Eq. (31). p*^: Slopes which give 
extrema of r{pr) (Eq. (40)). 



Fig. 21 is a schematic diagram of the r-pr curve, where Ces{4>) is assumed to be 
negative, r is given by r = dfcs{p)/dpr, and Pr is the surface slope along the line PP' 
(Fig. 7). Let us define p*_^ so that dr{pr)/dpr = 0. We then have 

1 



-3Seff(0) ± v/95eff(0)2 - 2AAM)CM) ■ (40) 

When yleg(0) > and i?eff(0) < 0, for < p < pi, a vicinal surface with a regular 
train of steps is in an equilibrium state, for pi < p < it is metastable, and for 
Pr+ < P < Po (Po = V^), it is thermodynamically unstable. 

For the temperature region Tf 2 < T < Tji, we have ?2g^''(7r/4) > in §3.3 and 
7^ > in §2.4, so that Acff(7r/4) > 0. From Eq. (36) and Eq. (37), we have 

" ' ' ' " 4A,ff(vr/4)' 16Aeff(vr/4)7i(7r/4)2- 

By comparing Eq. (41) with the PWFRG results in §2.4, we see that the non- 
universal shape exponents for the p-RSOS model agree with the values obtained using 
Eq. (41). Moreover, by comparing the amplitudes of Eq. (41) with the PWFRG results 
in §2.4, we obtain /3Aeff(7r/4) = (6.22 ± 0.06) x 10-^ and /37i(7r/4) = 1.39 ± 0.03. 
This value of /37i(7r/4) agrees well with /97(7r/4)ising (Eq. (18)). Consequently, in the 
Pr ^ limit, step droplets "evaporate" and dissociate into individual steps with n = 1. 
Substituting the value of PAesi'^'/'i) and the value of tiq^-* = 2.26 obtained from Fig. 19 
into Eq. (32), we have /3j[^\7i/4:) = (2.75 ±0.05) x lO'^. The small value oi (3j[^\Tr/4:) 
is consistent with experimental observations on Si(113)[93]. 
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5.2. T> Tf^i 

Prom the results in §2.4, a vicinal surface with small Pr shows the GMPT universal 
behavior of single steps. Near pr ^ 0, therefore, Aes(n/4:) — 0, and then 7^^^(7r/4) is 
considered to be zero. Hence, we have from Eq. (37), Eq. (36) and Eq. (31), 

On = 3/2, Ot - 26^ = 3, (42) 

By comparing Eq. (43) with the curve calculated using the PWFRG method, we have 
/3Seff(7r/4) = 0.281±0.008. /3Si(7r/4) at feT/e = 0.37 is estimated as /35i(7r/4) = 1.206 
by use of the step stiffness of the Ising model (Eq. (18)). Then, from Eq. (33) with 
= 1.45 (Pig. 19) and with 7? V/4) = 0, we have ^7^ V/4) = -0.88 ± 0.08. 



5.3. T< Tf^2 

In the low temperature region, the first-order shape transition occurs at the (001) facet 
edge on the ECS profile. If we consider Aefr(7r/4) < and Besin/A) < (Pig. 21), the 
first-order shape transition can be understood in terms of the non-GMPT vicinal surface 
free energy fes{p) (Eq. (31)). Prom the EPS results (§2.3), we have 7j"^''(7r/4) < 0, and 
from the Monte Carlo calculation of (n) at k-^T /e — 0.35, we have u^q'' > 0. These 
results are consistent with 74eff(7r/4) < 0. 



6. Summary and discussion 

In §2, we present statistical mechanical calculations on the vicinal surface of the 
restricted solid-on-solid (RSOS) model with a point-contact type step-step attraction 
(p-RSOS model) (Pig. 2). Applying the product wave-function renormalization group 
(PWPRG) method, which is a variant of the density matrix renormalization group 
(DMRG) method, to the transfer matrix for the p-RSOS model, we calculate the reduced 
equilibrium crystal shape (ECS) of the p-RSOS model (Pig. 6) and the surface gradient 
P = {Px,Py) as a function of X (= (3rix = —(3Xx, (3 = 1/k-QT) and Y (= (3riy = —/3\y) 
(Pig. 9). 

We obtain the first-order shape transition around the (111) facet on the ECS profile 
for T < Tf^i. Por T < T/_2, a first-order shape transition at the (001) facet edge is also 
observed on the ECS profile, where the (001) facet directly contacts the (111) facet 
(Pig. 10). By analyzing the PWPRG results, we obtain the equilibrium facet shape 
(Pig. 10) and the step tension; for the vicinal surface tilted towards the (110) direction, 
we obtain the step stiffness, the non-GMPT shape exponents (Pig. 11) and the non- 
GMPT amplitudes (Pig. 12). 
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In §3, in order to elucidate the origin of the non-GMPT behavior, we study step 
droplets formed by thermal step bunching on the vicinal surface tilted towards the (110) 
direction for the p-RSOS model. We calculate the vicinal surface free energy from the 
Andreev surface free energy obtained from the PWFRG calculations (Fig. 14). Step 
bunching near equilibrium caused by a singularity in the surface free energy (Fig. 15) is 
demonstrated using Monte Carlo simulations with a simple Metropolis algorithm (Fig. 
16, 17). To obtain a clear image of local step droplets, we demonstrate giant step 
formation similar to step faceting for T < From a microscopic point of view, 

these giant steps are formed from sticky steps due to the step-step attraction. For 
Tf,2 ^ T < T/ i, the giant steps partially dissociate due to entropic repulsion. Using a 
long-duration Monte Carlo simulation, we calculate the slope dependence of the mean 
step number (n) in a giant step (Fig. 19), and find that (n) increases linearly with 

In §4, we derive a |p|-expanded form of the non-GMPT vicinal surface free energy 
fes{p) (Eq- (31))- In Eq. (31), the |pp term appears. In the derivation, the concept 
of step droplets (boson n-mers. Fig. 20) and knowledge of the slope dependence of (n) 
obtained using the long-duration Monte Carlo calculation in §3.3 are crucial. 

In §5, the results in §2 are consistently reproduced by thermodynamical calculations 
based on the non-GMPT vicinal surface free energy fes{p)- By comparing the results 
in §2 with the results obtained using fcs{p), we obtain information on the derivative 
of 7„. d{'~fn/n)/dn\n=i is small and positive in the temperature range < 7" < T/^i, 
whereas it is negative for T < Tf 2 and zero for T > Tji. d'^{'-)n/n)/dn^\n=i < for 
T pa Tf^i and the step droplets arc dispersed in the limit \p\ — )> 0. 

The first-order transition on the ECS profile has been studied theoretically by 
several authors[26, 27]. Rottman and Wortis[26] calculated the interface tension using 
a 3D cubic Ising model with both nn and next nearest neighbor (nnn) interactions 
between spins by means of the mean field approximation. They found the first- 
order shape transition around the (001) facet at low temperature for the negative nnn 
interactions. However, they did not discuss the shape exponent for the temperature 
range < T < Tt (using their terminology). Jayaprakash et a/. [27] studied the 
vicinal surface of the interacting terrace-step-kink (TSK) model with long-range step- 
step attractions corresponding to attractive dipolar interactions. They determined the 
surface free energy using mean field calculations, and showed that the surface free energy 
f{p) has the form of the GMPT (Eq. (1)), and that the long-range attractions changes 
B. When B becomes negative, the facet edge causes the first-order shape transition. 
This explanation for the first-order shape transition is limited to the case of attractive 
long-range step-step interactions, because in the case of short-range attractions, B does 
not change. Hence, the explanation based on step droplets described in §4.2 is required. 

In real systems, there are a variety of different surface effects that occur, as described 
in §3.3.1. Of these, elastic interactions among steps are the most important source of 
long-range step-step repulsion. The TSK model with a long-range step-step repulsion 
of the order of ~ 1//^, where I is the inter-step distance, is known to have the GMPT 
universal free energy[27, 31] in a homogeneous system, and the system has the GMPT 
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universal shape exponents on the ECS. The work of Shenoy et a/. [81] suggests that step 
bunching occurs as m-mers of particles in the TSK model with a long-range step-step 
repulsion of the order of ~ 1//^ and a short-range step-step attraction. In this case, if 
m depends on p, then the non-GMPT shape exponents on the ECS are expected from 
Eq. (31)-Eq. (34). 

Paulin et a/. [94, 53] considered a surface system with an elastic interaction of 
~ 1/Plnl in order to study the dynamical behavior of the single-step to double-step 
transition on the Si(lOO) surface. Such an elastic interaction is thought to change 
the absolute values of the step tension, stiffness, and interaction coefficient (using 
our terminology), as in the case of ~ 1//^ [31]. As long as a homogeneous system 
is realized ("in phase meandering" using their terminology), an elastic interaction of 
~ 1//^ In / probably does not change the shape exponents. However, in the case of "out 
of phase meandering", since inhomogeneity is suggested, non-GMPT shape exponents 
are expected due to the step droplets as explained in §4.2. 

Recently, a first-order like shape change for ^He around the (0001) facet has been 
reported [95]. The strong anisotropy in the surface stiffness measured in the experiment 
suggests the existence of a singularity in the surface free energy. The results of the 
present study and those in our previous reports[54, 55] are consistent with this finding. 

7. Conclusion 

In the present study, we investigate a vicinal surface tilted towards the (110) direction 
near a (001) facet on the ECS of the p-RSOS model. For T/,2 < T < T/,i, the 
shape exponents on the ECS have non-GMPT values such as ^„ = 2 for the normal 
direction and 9t — 4: for the tangential direction. The origin of the non-GMPT shape 
exponents is the formation of step droplets ( "giant steps" ) with different sizes. The non- 
GMPT expression for the vicinal surface free energy fcs{p) is derived. In the derivation, 
knowledge of the \p\ dependence of the mean step-droplet size is crucial. The results 
obtained by statistical mechanical calculations using the p-RSOS model are successfully 
reproduced using thermodynamical calculations based on fes{p)- 
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Appendix A. Exact expressions for interface tension and interface stiffness 
on a two-dimensional square Ising model 

The exact expression for the 2D ECS for a nn square Ising model[68]-[71] is written 
using the angle (f) in Fig. 7(b) as follows: 

D{X,,Y,) = 0, 



dD{X,Y) 
dY 



{x,Y)={Xc,Yc) = tan( 



dD{X,Y) 
dX ' 

cosh^(/9e 



(X,y) = (Xe,Yc)> 



where D{X,Y) is 

D{X, Y) = cosh(X) + cosh(r) 

sinh(pej 

Hence, the equation for the ECS on the 2D square nn Ising model becomes 

cosh^(/3e) 



(A.l) 



(A.2) 



cosh(Xc) + cosh(yc) = 



sinh(/3e) ■ ^^-^^ 
Let us denote the anisotropic interface (step) tension for an interface with a mean 
running direction parallel to the hne QQ' by 7(0) (Fig. 7), and let us consider the 
interface stiffness 7(0) = 7(0) + 9^7(0)/(?0^. The exact expression for the interface 
stiffness for the 2D nn Ising model is derived by use of D{X,Y) (Eq. (A.l)) with Eq. 

(13), dL'(X(0),y(0))/d0|(x,y)=(x.,y.) = and d'D{X{(l>),Y{(l>))/d(l>\x,Y)=ix.,Y.) = 
as follows: 



Ising 



dD{X,Y) , , dD{X,Y) . ■ 
cos H — sm 



dX 



dY 



X 



-2 



d'D{X,Y) . , d'D{X,Y) 
-qJ^2 sin <^ + cos 

d^D{X,Y) 
dYdX 



QY2 

-\ -1 



sin cos ( 



(A.4) 



(x,y)=(Xe,yc) 

Appendix B. Derivation of the non-GMPT vicinal surface free energy 



We investigate the n dependence of fn{Pn) for Pn ~ 0. Pn is expressed as p„ = Ng^n/ L — 
\p\/dn, where A^s,n is the total number of giant steps with size n, L is the length of the 
projected area of the vicinal surface, and dn is the height of the giant step. L is assumed 
to be sufficiently large so that we may take the thermodynamic limit. Substituting the 
expressions for and dn = ndi into Eq. (27), we have 



Up) = /(o) + 



7n(0), , , Bni<p) 



din 



\P\ + 



dfn^ 



P\ + 



d\n^ 



p\' + 0{p'>). 



(B.l) 



lri\ 



is considered to be approximately equal to n7i(0), we expand 7„(0)/n with 



Since 

respect to n around n = 1 as follows: 



71 (0) + lY'mn - 1) + ^7^(0) (n - 1)^ + • 



n 



(B.2) 



n=l 
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(a) (b) 

Figure CI. (a) Choice of X- and F-axes at P on the facet contour, (b) Q lies near P 
on the facet contour. The "normal distance" AX' is related to "tangential distance" 
Ay as AX' = (Ay)2/2i?, where R is the radius of curvature at P. 



Assuming the universal relation Eq. (17), where 7„(0) represents the stiffness of the 
giant step 7n(0) = 7n(0) + ^^7n(0)/50^, we expand %{(f)) instead of Bn{(f)). That is, 



7n(0) 



71 (0) + l['\4>){n - 1) + U?\cP){n + - 



n ' ' ' ' ' ' 2 



7i Wl = 



(B.3) 



9n™ 

Substituting Eq. (B.2) and Eq. (B.3) into Eq. (B.l), we have Eq. (28). 

Appendix C. Thermodynamic derivation of the shape exponents 

Eq. (36) - Eq. (39) are derived from the extended vicinal surface free energy jg{p) 
(Eq. (35)). In the following manner, we consider the case where the first-order shape 
transition docs not occur on the contour of the (001) facet. 

From the thermodynamics of the ECS[8], the coordinates {X^Y^ on the reduced 
ECS are obtained from /(p) as follows[28]: 

x = /jMz^, y = /?^, (0.1) 

where (3 = l/k-oT. Then the reduced coordinates {X.Y) on the ECS are written 

X = X,(0) + IpI^-^ [/3C5c(0) cos</> ~ PB[{^) sin0] , 
Y = ye(0) + IpI^-' [KBci^) sin0 + /3i?[(0) COS0] , 

B[{ct>) = dB^{ct>)ldct>, (C.2) 

where Xc and Yc are given as follows: 

Xc = /3'y{(f)) COS0 — /3j'{(f)) sin0, 
Yc = /37(0) sin0 + /37'(0) COS0, 

7'(0) = a7(0)/a0. (C.3) 



Non-universal equilibrium crystal shape 



28 



Let us choose the X- and Y- axes so that the F-axis is paraUel to the tangential 
hne of the facet contour at P (Fig. CI). With this choice of coordinate system, wc have 
(/) = at P and 7'(0) = (Eq. (C.3)). Along the tangential line {X = Xc{0) = 0), 
and |p| are not independent but are constrained to satisfy 
1 



-7(0)0^ + CSc(0)|p|^-^ = 0, 



(C.4) 



which is derived by expanding (C.2) and (C.3) with respect to and |p| (|p| << 1 and 
\(f)\ « 1, near the point P). Combining (C.4) with (C.2) and (C.3), we obtain 







^(oT' 



_(^bT)^AF^ 
2C7(0)i?c(0) 



c-i 



along the tangential line PQ in Fig. CI, where AY 
f{Xc{0), AY + Yc{0))/kBT is expanded to give 



(C.5) 

Y - Yc. Near P, Z = 



AZ^ 



[7(0) + '^YjoW 

kT,T 



\P\ + 



P\ 



Ar|p|0-x,(o)|p|(i--02), 



(C.6) 



where AZ — Z{X, Y) — Z{Xc, Yc). Substituting (C.5) into (C.6), we obtain expressions 
for the tangential shape exponent 9t and the tangential amplitude ^t(0) as follows: 

Az = -A(o)|Ay|^*, 



A(o) 



1 



{kBT) 



C+i 



c-i 



^, = 2C/(C-i) 



(C.7) 



.2CSc(0)7(0)C. 

Similarly, along the normal line {Y — Yc{0) = 0), we obtain expressions for the normal 
amplitude Ani4>) and the normal shape exponent 9n as follows: 

AZ = -AnmAXf", (Xe(0) < X) 



An{0) 



On 



kBT 

LC5c(0)J 



C-1 



^n = C/(C-l)- 



From Eq. (36) and Eq. (37), we obtain a scaling relation like 
Ot — 2^„ 

when 7(0) 7^ 0. From Eq. (37) and Eq. (36), we obtain the relation 



AX 



MO) 
A(o) 



AY' = MLay\ 
27(0) 



(C.8) 



(C.9) 



(C.IO) 



This expression is also obtained geometrically from Fig. Cl(b) using the radius of 
curvature R such that AX' = AY'/{2R) = A;BTAy 7(27(0)) [66] and AX' = AX for 
small AX' and AY. 
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